Overview

I know I billed this as “Ditching ArcGIS” in the poll, but I 1) am not an ArcGIS expert (despite using it for many years) and 2) there may still be things that you must do in ArcGIS to get your work done. If that’s the case, there is the R - ArcGIS Bridge that makes it possible for these two analysis tools to play with one another, BUT I won’t discuss that here…

Configure chunk options

# Set global knitr chunk options
if (.Platform$OS.type == "unix") {
  # Do not specify Cairo device for MacOS
  knitr::opts_chunk$set(warning = F, message = F, out.width = '100%', fig.align = 'center')
} else {
  knitr::opts_chunk$set(warning = F, message = F, out.width = '100%', fig.align = 'center', 
                        dev.args = list(type = "cairo"))
}

Load libraries

You may need to install some packages (and their dependencies) that you don’t already have. The chunk below will do that for you automatically using the load_pkgs function (hopefully you don’t mind).

# List packages required to run the script -------------------------------------
pkgs <- c("tidyverse","ggmap","sf","cowplot","here","devtools",
          "knitr","ggrepel","rgeos","FedData","raster","mapview",
          "rnaturalearth","rnaturalearthdata","shadowtext","leaflet",
          "leaflet.extras","htmltools")
# Install and load all CRAN packages provided from a character vector
load_pkgs = function(pkgs) {
  new_pkgs = pkgs[!(pkgs %in% installed.packages()[ ,'Package'])]
  if (length(new_pkgs) > 0) install.packages(new_pkgs,repos = "http://cran.cnr.berkeley.edu/")
  invisible(lapply(pkgs,function(x)
    suppressPackageStartupMessages(library(x,character.only = T))))
}
# Load packages
load_pkgs(pkgs)

# Load dev packages -------------------------------------------------------
# Install and load rnaturalearthdata package from github
if ("rnaturalearthdata" %in% installed.packages()[ ,'Package'] == F) {
  devtools::install_github("ropenscilabs/rnaturalearthdata")
}
library(rnaturalearthdata)

# Install and load rnaturalearthhires package from ropensci
if ("rnaturalearthhires" %in% installed.packages()[ ,'Package'] == F) {
install.packages("rnaturalearthhires",
                 repos = "http://packages.ropensci.org",
                 type = "source")
}
library(rnaturalearthhires)

# Create output directories if missing
dir.create(here("Figs"))
dir.create(here("Output"))

A VERY basic map

Map with base graphics

Let’s create a very basic “map”, showing the location of acoustic samples along the West Coast. The first uses base graphics (yuck!) and a cartesian coordinate system.

# Load some ship track data
cps <- read_csv(here("Data/cps_nav.csv"))

# A basic 
plot(cps$long, cps$lat)

# plot(cps$long, cps$lat, cex = cps$cps.nasc)
# plot(cps$long, cps$lat, cex = sqrt(cps$cps.nasc/1000))

Map using ggplot2

Here are the same data plotted using ggplot2.

# A very basic ggplot2 scatter plot of the same data
ggplot(cps, aes(long, lat)) + geom_point()

A slightly less basic map

Let’s set a few parameters to make the map more spatially accurate. But it’s a long way from a useful “map”…

# Add some more aesthetics and coordinate controls, maybe some labels

# Read locations
locations <- read_csv(here("Data/locations.csv")) %>% 
  filter(group %in% c("city", "landmark"))

# Refine the "map"
ggplot(cps, aes(long, lat)) + 
  # geom_point() +
  # geom_path(aes(group = transect, colour = factor(transect))) +
  geom_point(aes(group = transect, colour = factor(transect), size = cps.nasc)) +
  geom_point(data = locations, aes(long, lat)) +
  # geom_text(data = locations, aes(long, lat, label = name)) +
  # coord_quickmap() +
  coord_map(projection = 'azequalarea') +
  # coord_map(projection = 'azequalarea', xlim = c(-128,-123), ylim = c(45,50)) +
  theme_bw()

Canned mapping options

Before we go farther, let’s set some map boundaries based on either our data (cps$lat and cps.long) or manually (mb.lat and mb.long).

# Define lat and long bounds for west coast map
wc.lat  <- range(cps$lat)  #c(32, 52)
wc.long <- range(cps$long) #c(-130, -116)

# Define lat and long bounds for Monterey Bay
mb.lat  <- c(36.5, 37)
mb.long <- c(-122.2, -121.7)

ggmap

ggmap is a mashup of ggplot and online mapping services. There are numerous map sources available to ggmap (e.g., Google, OpenStreetMaps, Stamen). We’ll look quickly at Stamen and Google maps. Withing each of those, there are multiple options (e.g., Google has satellite, terrain, etc.). Each has different arguments to specify the map boundaries.

Stamen maps

Create a basic Stamen map (Toner option). Look, a map (with no data)!

# Set west coast boundaries for stamen maps
wc.bounds.stamen <- c(left = min(wc.long), bottom = min(wc.lat),
                      right = max(wc.long), top = max(wc.lat))
# Download stamen map of west coast; zoom = 6 seems good
wc.map.stamen.toner <- get_stamenmap(wc.bounds.stamen, zoom = 6, maptype = "toner-lite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw()
# Display the stamen map
wc.map.stamen.toner

Add our data to the map. Look, a map (with data)!

# Add layers to map
wc.stamen1 <- wc.map.stamen.toner + 
  geom_point(data = cps, aes(long, lat, group = transect, colour = factor(transect), size = cps.nasc),
             show.legend = F) +
  geom_text(data = locations, aes(long, lat, label = name)) +
  ggtitle("Basic options")

wc.stamen2 <- wc.map.stamen.toner + 
  geom_point(data = cps, aes(long, lat, group = transect, colour = factor(transect), size = cps.nasc),
             show.legend = F) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text(data = locations, aes(long, lat, label = name), 
            colour = "red", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) + 
  ggtitle("Formatted labels")

wc.stamen3 <- wc.map.stamen.toner + 
  geom_point(data = cps, aes(long, lat, group = transect, colour = factor(transect), size = cps.nasc),
             show.legend = F) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text_repel(data = locations, aes(long, lat, label = name), 
                  size = 2, segment.colour = "black", segment.alpha = 0.5) +
  ggtitle("Formatted and repelled lables")

# Combine maps
wc.grid.stamen <- plot_grid(wc.stamen1, wc.stamen2, wc.stamen3, nrow = 1)

# Save map
ggsave(wc.grid.stamen, filename = here("Figs/wc_map_stamen_toner.png"), width = 10, height = 7)
# Print map
include_graphics(here("Figs/wc_map_stamen_toner.png"))

Google maps

# Reduce label list
label.list <- c("Monterey Bay","San Francisco","Cape Flattery","Crescent City",
                "Newport","Point Conception","Cape Mendocino","Columbia River",
                "Cape Blanco","Bodega Bay","Westport","Fort Bragg",
                "Morro Bay","Long Beach","Cape Scott","San Diego")

locations <- filter(locations, name %in% label.list)

# Set west coast boundaries for stamen maps
wc.bounds.google <- c(mean(wc.long),mean(wc.lat))
# Download stamen map of west coast
wc.map.google <- get_map(location = wc.bounds.google, zoom = 4, source = "google", maptype = "satellite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
  xlim(wc.long) + ylim(wc.lat)

wc.google1 <- wc.map.google +
  geom_point(data = cps, aes(long, lat, group = transect, colour = cps.nasc, size = cps.nasc)) +
  scale_colour_gradientn(colors = rev(rainbow(10))) +
  geom_point(data = locations, aes(long, lat), colour = "white") +
  geom_shadowtext(data = locations, aes(long, lat, label = name), 
                  colour = "white", bg.color = "black", size = 3, 
                  nudge_x = 0.25, hjust = 0, angle = 45, fontface = "bold.italic") +
  ggtitle("Shadow text, scale color & size")

wc.google2 <- wc.map.google +
  geom_point(data = cps, aes(long, lat), colour = "gray50", size = 0.25, alpha = 0.75) +
  geom_point(data = filter(cps, cps.nasc > 0), 
             aes(long, lat, group = transect, colour = cps.nasc, size = cps.nasc), alpha = 0.75) +
  scale_colour_viridis_c(option = "magma") +
  geom_point(data = locations, aes(long, lat), colour = "white") +
  # Configure legend guides
  guides(colour = guide_legend(), size = guide_legend()) +
  geom_shadowtext(data = locations, aes(long, lat, label = name), 
                  colour = "white", bg.color = "black", size = 2, 
                  nudge_x = 0.25, hjust = 0, angle = 45, fontface = "bold.italic") +
  ggtitle("Viridis colours, grouped legend")

# Arrange plots
wc.grid.google <- plot_grid(wc.google1, wc.google2, nrow = 1)

# Save map
ggsave(wc.grid.google, filename = here("Figs/wc_map_google.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_google.png"))

# Set west coast boundaries for stamen maps
mb.bounds.google <- c(mean(mb.long),mean(mb.lat))
# Download stamen map of west coast
mb.map.google <- get_map(location = mb.bounds.google, zoom = 10, 
                         source = "google", maptype = "satellite") %>% 
  ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
  xlim(mb.long) + ylim(mb.lat)

mb.google1 <- mb.map.google +
  geom_point(data = filter(locations, name != "Monterey Bay"), aes(long, lat), colour = "white") +
  geom_text(data = filter(locations, name != "Monterey Bay"), aes(long, lat, label = name), 
            colour = "white", size = 4, vjust = 0, nudge_y = 0.01)
# Save map
ggsave(filename = here("Figs/mb_map_google.png"), width = 7, height = 7)
# Print map
include_graphics(here("Figs/mb_map_google.png"))

“Roll-your-own” mapping options

rnaturalearth

rnaturalearth for shoreline and country files. Below is a comparison of the medium- and large-scale country datasets.

# Import NaturalEarth coastlines and countries
coast_sf        <- ne_coastline(scale = "medium", returnclass = "sf")
coast_sf_lg     <- ne_coastline(scale = "large", returnclass = "sf")
countries_sf    <- ne_countries(scale = "medium", returnclass = "sf")
countries_sf_lg <- ne_countries(scale = "large", returnclass = "sf")

# Display countries
plot(countries_sf)

# View unique regions
sort(unique(countries_sf$region_wb))
## [1] "Antarctica"                 "East Asia & Pacific"       
## [3] "Europe & Central Asia"      "Latin America & Caribbean" 
## [5] "Middle East & North Africa" "North America"             
## [7] "South Asia"                 "Sub-Saharan Africa"
sort(unique(countries_sf$subregion))
##  [1] "Antarctica"                "Australia and New Zealand"
##  [3] "Caribbean"                 "Central America"          
##  [5] "Central Asia"              "Eastern Africa"           
##  [7] "Eastern Asia"              "Eastern Europe"           
##  [9] "Melanesia"                 "Micronesia"               
## [11] "Middle Africa"             "Northern Africa"          
## [13] "Northern America"          "Northern Europe"          
## [15] "Polynesia"                 "Seven seas (open ocean)"  
## [17] "South-Eastern Asia"        "South America"            
## [19] "Southern Africa"           "Southern Asia"            
## [21] "Southern Europe"           "Western Africa"           
## [23] "Western Asia"              "Western Europe"
# Select only certain regions to speed plotting
na_sf    <- filter(countries_sf, subregion %in% c("Northern America","Central America")) 
na_sf_lg <- filter(countries_sf_lg, subregion %in% c("Northern America","Central America")) 

# Create West Coast map (medium scale)
wc.ne <- ggplot() + 
  geom_sf(data = na_sf) +
  geom_point(data = locations, aes(long, lat)) +
  geom_text(data = locations, aes(long, lat, label = name), size = 4) +
  scale_x_continuous(name = "Longitude", limits = wc.long) + 
  scale_y_continuous(name = "Latitude", limits = wc.lat) +
  theme_bw() + 
  coord_sf(xlim = c(-123.5,-121.5), ylim = c(36.5,38.5)) +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5),
          legend.position =  c(0,0),
          legend.justification = c(0,0),
          legend.background = element_blank(),
          legend.key = element_blank(),
          panel.background = element_rect(fill = alpha("lightblue", 0.5)),
          plot.title = element_text(hjust = 0.5),
          panel.grid.major = element_line(color = "white")) +
  ggtitle("Medium scale") 

# Create West Coast map (large scale); and with shadowtext labels
wc.ne.lg <- ggplot() + 
  geom_sf(data = na_sf_lg) +
  geom_point(data = locations, aes(long, lat)) +
  geom_shadowtext(data = locations, aes(long, lat, label = name), colour = "white",
                  size = 4, bg.colour = "black") +
  scale_x_continuous(name = "Longitude", limits = wc.long) + 
  scale_y_continuous(name = "Latitude", limits = wc.lat) +
  theme_bw() + 
  coord_sf(xlim = c(-123.5,-121.5), ylim = c(36.5,38.5)) +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5),
          legend.position =  c(0,0),
          legend.justification = c(0,0),
          legend.background = element_blank(),
          legend.key = element_blank(),
          panel.background = element_rect(fill = alpha("lightblue", 0.5)),
          plot.title = element_text(hjust = 0.5),
          panel.grid.major = element_line(color = "white")) +
  ggtitle("Large scale")

# Create map grid
wc.grid.ne <- plot_grid(wc.ne, wc.ne.lg, nrow = 1)

# Save map
ggsave(wc.grid.ne, filename = here("Figs/wc_map_natEarth.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_natEarth.png"))

Interactive maps

Leaflet

Leaflet is an open-source JavaScript library for creating interactive maps. The Leaflet for R site is indispensible for creating these maps. We’ll show mapview below once we’ve converted the data to a spatial object (sf in this case).

# Select tile to use for Leaflet map
# Some good options include CartoDB.Positron, Stamen.Terrain, Esri.WorldImagery
leaflet.tile <- "Stamen.Terrain" 

# A very basic Leaflet map
leaflet() %>% 
  # addTiles() %>% 
  addProviderTiles(leaflet.tile) %>% 
  addMarkers(data = locations, ~long, ~lat, popup = ~as.character(name), label = ~as.character(name))

GIS tasks using the sf package

I recently (and fortuitously) had a real-life GIS problem that I needed to solve. Every three years, we must renew our CDFW Scientific Collection Permit, and to do so, we must identify how many trawls we conducted in CA State waters, and how many we did in MPAs (if any). Below are the steps I took to automate this process (data processing steps not shown).

Reading shapefiles

I have some shapefiles of the CA state waters and CA marine protected areas (MPAs). I need to read them for plotting later. You often must know the datum and/or projection. These can be defined using proj4 (I’m not going to talk about this) or the European Petroleum Survey Group (EPSG) or coordinate reference system (CRS) codes. You can search for these codes at the Spatial Reference List. Some commonly used ones are 4326 for WGS84 (for GPS data) and NAD83 (4269, for natural earth). You may want more spatially conservative projections for small scales or where preserving spatial accuracy is critical.

# Get map data
  na_sf <- ne_countries(scale = "large", returnclass = "sf") %>%
    filter(name %in% c("United States", "Mexico", "Canada")) %>% 
    st_transform(4269) # NAD83
  
  # Get CA State Waters shapefile
  ca_waters <- st_read(here("GIS/MAN_CA_StateWater.shp")) %>% 
    st_transform(4326) # WGS84
## Reading layer `MAN_CA_StateWater' from data source `D:\CODE\R\R-users\rGIS\GIS\MAN_CA_StateWater.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -384564.5 ymin: -605050.3 xmax: 270524.5 ymax: 450642.2
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs
  # CA MPAs
  ca_mpas <- st_read(here("GIS/MPA_CA_Existing.shp")) %>% 
    st_transform(4326) %>% 
    mutate(MPA = paste(NAME, Type))
## Reading layer `MPA_CA_Existing' from data source `D:\CODE\R\R-users\rGIS\GIS\MPA_CA_Existing.shp' using driver `ESRI Shapefile'
## Simple feature collection with 121 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -373440.2 ymin: -590425.8 xmax: 259220.6 ymax: 259522.3
## epsg (SRID):    NA
## proj4string:    +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs

Converting data frames to simple features

We want (need?) to load the haul locations and convert to simple feature before doing GIS stuff.

load(here("Data/haul_df.Rdata"))

# Convert to sf; CRS =4326 (WGS84)
haul_sf <- st_as_sf(haul, coords = c("long","lat"), crs = 4326) %>% 
  mutate(long = map_dbl(geometry, ~st_centroid(.x)[[1]]),
         lat  = map_dbl(geometry, ~st_centroid(.x)[[2]]))

# We'll also convert our locations to sf
locations <- locations %>% 
  st_as_sf(coords = c("long","lat"), crs = 4269) %>% 
  mutate(
    long = map_dbl(geometry, ~st_centroid(.x)[[1]]),
    lat  = map_dbl(geometry, ~st_centroid(.x)[[2]]))

Create a static base map

# Set padding around data  
  bbox.loc <- st_bbox(haul_sf)
  
  # Create base map
  base.map <- ggplot() +
    # Plot high-res land polygons
    geom_sf(data = na_sf, fill = "white", color = "tan4") +
    # Plot landmarks
    geom_sf(data = locations, size = 2, colour = 'tan4') +
    geom_shadowtext(data  = locations, aes(long, lat, label = name), colour = 'gray20', size = 2,
              fontface = 'bold', hjust = 0, nudge_x = 0.2, nudge_y = 0.05, 
              angle = 25, bg.color = "white") +
    # Format axes and titles
    xlab("Longitude") + ylab("Latitude") + 
    coord_sf(
      xlim = c(bbox.loc$xmin,bbox.loc$xmax),
      ylim = c(bbox.loc$ymin,bbox.loc$ymax)) +
    theme_bw() + 
    theme(axis.text.y = element_text(angle = 90, hjust = 0.5),
          legend.position =  c(0,0),
          legend.justification = c(0,0),
          legend.background = element_blank(),
          legend.key = element_blank(),
          panel.background = element_rect(fill = alpha("lightblue", 0.5)),
          plot.title = element_text(hjust = 0.5),
          panel.grid.major = element_line(color = "white"))

Add some data to the base map

Let’s add our CA state waters and MPAs to the map, along with our trawl haul data.

# Add CA state waters layer
haul.map <- base.map + 
  geom_sf(data = ca_waters, colour = "red", fill = NA) +
  geom_sf(data = ca_mpas, aes(fill = Type), colour = "gray50") +
  geom_sf(data = haul_sf, colour = "blue", shape = 21) +
  coord_sf(
      xlim = c(bbox.loc$xmin,bbox.loc$xmax),
      ylim = c(bbox.loc$ymin,bbox.loc$ymax))

# Show map
haul.map

Find an intersection

# Find hauls in CA waters
haul.ca <- st_intersection(haul_sf, ca_waters)

haul.mpa <- st_intersection(haul_sf, ca_mpas) %>% 
  mutate(key = paste(cruise, ship, haul, collection))

Intersection results; static map

Add the results of the intersection to the map.

haul.map + 
  geom_sf(data = haul.ca, colour = "green") +
  geom_sf(data = haul.mpa, colour = "yellow") +
    coord_sf(
      xlim = c(bbox.loc$xmin,bbox.loc$xmax),
      ylim = c(bbox.loc$ymin,bbox.loc$ymax))

Intersection results: Leaflet

# Configure palette for MPAs
factpal <- colorFactor(topo.colors(10), ca_mpas$MPA)

haul.out <- filter(haul, !key %in% haul.ca$key)

# Create leaflet map
leaflet() %>% 
  addProviderTiles(leaflet.tile) %>% 
  addPolygons(data = ca_waters, weight = 2, fillColor = "transparent") %>% 
  addPolygons(data = ca_mpas, color = "gray50", weight = 2, fillColor =  ~factpal(MPA), fillOpacity = 1,
              label = ~htmlEscape(MPA)) %>% 
  addCircleMarkers(data = haul.out, radius = 3, color = "gray50", stroke = FALSE, fillOpacity = 0.75,
                   label = ~htmlEscape(key)) %>% 
  addCircleMarkers(data = haul.ca,  radius = 5, color = "red", stroke = FALSE, fillOpacity = 0.75,
                   label = ~htmlEscape(key)) %>% 
  addMarkers(data = haul.mpa, popup = ~key, label = ~htmlEscape(key))

Intersection results: Mapview

mapview(ca_waters) +
  mapview(ca_mpas) +
  mapview(haul_sf)

Data packages

tigris

Access to the TIGER/U.S. Census line data files on github.

elevatr

The elevatr package provides easy access to elevation data from a variety of sources. Not fully functional at the moment…

FedData

I’ve not had much time to explore this one.

FedData Github repo and tutorial with examples of various data sources (e.g., .

FedData Tutorial

Presently not executed; still under development

# Get a map of Moss Landing
ML.bbox <- get_map("Moss Landing, CA", zoom = 14) 
# Plot the map
ggmap(ML.bbox) 
# Get the bounding box from its attributes
bb <- attr(ML.bbox, "bb")
# Create an extent polygon from the Monterey Bay bounding box
ML.extent <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon,
                                            bb$ll.lat, bb$ur.lat), 
                             proj4string="+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
# Map the bounding box to check the results
ggmap(ML.bbox) + geom_polygon(aes(x=long, y=lat), size=3, color='purple',
                             data=ML.extent, fill=NA)
# Download National Elevation Data
ned_ML <- get_ned(template = ML.extent, label="ned_ML", 
                    res="1", force.redo = F, raw.dir = here("Data"))
# Convert raster to data frame for plotting in ggplot
ned_df <- rasterToPoints(ned_ML) %>% 
  data.frame() %>% 
  rename(lat = y,
         long = x,
         elev = grdn37w122_1)

# Map elevation
ggplot(ned_df, aes(long,lat,fill = elev)) + geom_raster() +
  geom_contour(aes(z = elev), alpha = 0.5, binwidth = 1) + 
  scale_fill_viridis_c(name = "Elevation (m)") +
  theme_bw()

Other potentially useful packages

More R-GIS resources

Simple Feaures for R: Blogs, presentations, vignettes for the sf package

A Tidy Approach to Spatial Data: Simple Features from our own Eric Anderson

GIS with R: an introduction by Francisco Rodriguez-Sanchez (Pakillo)

An Intro to GIS with R by Jess Sadler (covers sp, sf, and rnaturalearth, among other things).

Geocomputation with R

r-sptatial Blog